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FOREWORD 


The Industrial Mathematics Society is a professional organiza- 
tion whose objective it is to extend the understanding and application 
of mathematics in industry. Founded early in 1949, it marked the 
first organized attempt of any group to foster a closer relationship 
between mathematics and industry. 


It promotes the cause of mathematics in industry through 
monthly meetings at which technical papers are presented, through 
lectures by nationally prominent scientists and engineers, through 
panel groups which concentrate in specialized fields of knowledge, 
and through publication of worthwhile papers in its annual volume. 


Membership is open to engineers, scientists, and mathematicians 
who are concerned with the effective use of mathematics in industry. 
Annual dues are $5.00 per year including a subscription to the an- 
nual volume, Jndustrial Mathematics. Application forms and infor- 
mation may be obtained by writing to the Secretary, Industrial 
Mathematics Society, 100 Farnsworth, Detroit 2, Michigan. 
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Vice President, Robert A. Roggenbuck, Ford Motor Co., Engineering 
Staff 
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Data Processing for Numerical Machining * 


By E. C. Johnson 
Research Laboratories Division 
Bendix Aviation Corporation, Detroit 


INTRODUCTION 


Numerically controlled manufacturing has as its objective the 
production of finished parts directly from basic design information 
with a minimum of human intervention. The over-all concept is 
illustrated pictorially in Figure 1 in terms of the present Bendix 
system. As indicated, the system consists of two principal groups 
of equipment: the tape-preparation equipment, shown on the left, 
and the machine tool and its directly connected control equipment on 
the right. 
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TAPE PREPARATION MACHINING 


This paper is concerned primarily with the data-processing 
aspects of the Bendix system —more specifically, those activities 
associated with the first category of equipment just mentioned. This 
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equipment is physically distinct from the machine control unit and 
usually located inan office-type environment away from the machine 
area. Its function is to produce a record—in this case a punched 
plastic tape — suitable for controlling the motions of the machine. 

Involved as a first step in the operation is the organization and 
entry into a computer of data defining the part and the way in which 
it is to be machined. From this information, a sequence of cutter- 
center positions is generated, taking into account the required cutter 
offset. The result is a punched tape containing detailed instructions 
in coded form which enable the machine control unit to produce the 
proper relative motion between cutter and workpiece. 


INFORMATION REQUIRED BY THE MACHINE CONTROL UNIT 


Linear Interpolation 


The problems involved in tape preparation are determined in 
large measure by the nature of the information required by the 
machine control unit. In this respect, one of the more significant 
characteristics of the Bendix system is the ability to perform an 
elementary curve-fitting or interpolating operation at the machine. 
For example, with points given as indicated in Figure 2, the control 
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Figure 2. Straight-Line Approximation to a Curve 
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DATA PROCESSING FOR NUMERICAL MACHINING 3 
unit can in effect generate the connecting straight lines with no 
additional information. 

This ability to interpolate has important implications at the 
machine in terms of operating features and versatility. From the 
data-handling point of view, it likewise has important implications 
in greatly reducing the amount of information that must be trans- 
mitted to the machine. Data need be supplied to it only frequently 
enough that the resulting straight-line motions generate the desired 
part with acceptable tolerance. This type of information is relatively 
easy to produce with the aid of automatic computing equipment. 


Control Tape Format 


The information required for each straight-line motion is pro- 
vided to the machine control unit in terms of the differences in 
position along each of the coordinate axes, Ax, Ay, and Az, as in- 
dicated in Figure 2. This information is grouped on the control tape 
in blocks, one block for each of the elementary straight-line motions, 
and punched in channels which run lengthwise along the tape. The 
channels are numbered in conventional teletype equipment as shown 
at the bottom in Figure 3. The unnumbered channel near the center 
contains sprocket holes for manipulating the tape mechanically or 
for defining a line when the tape is read photoelectrically. 

In channel 1 is placed information which specifies A x, the mo- 
tion along the x axis. Channel 2 contains A y, and channel 3, Az, in 
Similar fashion. The last line of the block (near the numbered end) 
is used to indicate the sign, or direction, of motion by a punch for 
minus and not for plus. The remainder of the holes in each channel 
specify the magnitude of the motion desired. This is done in terms 
of basic units of machine motion, typically 0.0002 inch, using 
straight binary code. 

The contents of channel 4 is called the feed-rate number. This 
number determines the rate at which command signals are delivered 
to the servo drives and thereby controls the cutter speed. Channel 
5 is not read by the machine control unit. It is punched by the Bendix 
computer in preparing the tape to make possible the reading back of 
tape into the computer for checking or for modification. 

Channel 6 defines the length of tape associated with each block 
by means of a punch at the end of the biock. Each biock, then, is 
only as long as necessary to contain the maximum number of binary 
digits required for any of the numbers Ax, Ay, Az, or feed rate. 
Channel 7 is used for checking purposes to help insure that the tape 
is punched and read correctly. Holes are punched or not punched in 
this channel in order to make the total number of holes in each line 
across the tape an odd number, except for channel 5 which is not 
read at the machine. 
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The block described so far is a typical block of dimensional 
control information. Another kind of block (shown at the bottom in 
Figure 3) provides for handling auxiliary or on-off machine func- 
tions. An auxiliary-function block is indicated by punching both the 
end-of-block channel and the feed-rate channel in the same line, a 
combination that never occurs normally. Channels 1, 2, and 3 then 
specify which auxiliary function is desired. The block can be 
made as long as needed to handle any desired number of auxiliary 
functions. 


COMPUTATIONAL PROBLEMS 


The basic scheme described in the preceding section provides 
for moving the slides of a machine in a perfectly general way. It 
remains Only to so instruct the machine that such motion produces 
the desired effect in terms of the resulting workpiece. 
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DATA PROCESSING FOR NUMERICAL MACHINING 

Two principal computational problems are involved. The first 
results from the fact that motion of the machine slides is in effect 
motion of the cutter axis or center. In general, however, it is the 
periphery of the cutter that produces the surface defined in the part 
design information. The deduction of a suitable cutter-center path 
form knowledge of the desired finished surface is referred to as the 
cutter -offset problem. 

The second problem is that of generating sufficient points on a 
curve that the linear interpolation performed by the machine control 
unit produces the desired accuracy and surface finish. These 
problems are illustrated in Figure 4. Some additional problems of 
a less fundamental nature are those of computing feed-rate number 
and conversion of the usual decimal information into binary code for 
the Bendix machine control unit. 


APPROXIMATION 
TOLERANCE 


CUTTER 





These problems taken together are of sufficient magnitude to 
make desirable the use of a computer for all but the simplest of 
parts. Control tapes can, of course, be produced by means of a 
manual punch from hand-computed results. For complex parts, 
however, the economic success of numerical control will probably 
depend on making effective use of automatic computers of some sort 
in data preparation. 


THE BENDIX TAPE PREPARATION SYSTEM 


The Bendix approach to this problem makes use of a small but 
efficient general-purpose digital computer, equipped with specially 
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C. JOHNSON 
developed programs or routines. Figure 1 illustrates the principal 
elements in the system and the sequence of operations involved in 
tape preparation. 

The data required as input to the computer are of two kinds; 
dimensional data defining the part geometry and information con- 
cerned with the metal-cutting aspects of the job —cutter size, feed 
rates, sequence of cuts, etc. As in other systems of this type, a 
process sheet or manuscript is used as the means for organizing 
information for input to the computer. The process sheet itself is 
merely a printed form which serves as a work sheet for the process 
planner. However, the act of filling out the sheet automatically 
compiles and arranges in proper sequence the information needed 
by the computer. 

The next step in the process is that of producing from the pro- 
cess sheet a tape which can be read by the computer. This tape, re- 
ferred to as the process tape, contains exactly the same information 
as the process sheet. It is produced by a device called a Flexo- 
writer — basically an electric typewriter with a tape reader and tape 
punch attached. The punching of the process tape is done as the 
original handwritten process sheet is copied on the typewriter. 

Accuracy of the transcription can be checked by comparing the 
typewritten copy of the process sheet with the handwritten version. 
However, a more complete and automatic check can be made through 
use of the verifier. When the original tape is finished, it is inserted 


in the motorized tape reader connected to the Flexowriter. The 
Same process sheet is manually typed again to produce a second | 
tape while the first tape is read in synchronism. As each typewriter | 


key is struck, if the associated code agrees with that punched in the 
first tape, the tape punch operates and punches the verified code in 
the second tape. [If the codes do not agree, the punch does not oper- 
ate and the keyboard locks. The operator can then determine where 
the error lies and punch the second tape accordingly. 

At this point all manual effort involved in tape preparation is 
completed. The process tape now contains in condensed and coded 
form all necessary data and is fed into the computer by way ofa 
tape reader. The computer interprets the data provided by the 
process tape, carries out all the necessary computations, arranges 
the results in the required format, and punches the control tape 
exactly as required for input to the machine control unit. 

Figure 5 shows in more detail the Flexowriter used for punching 
the process tape from the handwritten process sheet. As evident in 
the Figure, the Flexowriter resembles an oversize typewriter with 
a tape punch and a tape reader mounted on the side. The tape 


reader provides a means for automatically producing a typewritten | 


copy of the contents of a process tape, if that should be desirable, 
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DATA PROCESSING FOR NUMERICAL MACHINING 


and it provides a means for duplicating both process tapes and con- 
trol tapes. In the rear and to the left is the tape reader used in the 
verification process. 





Figure 5, Photograph of Flexowriter-Verifier 


The computer used in the Bendix system, as already mentioned, 
is a small general-purpose digital computer, the Bendix G-15D, 
Figure 6. For the numerical machining application, a special input- 
output unit, the AN-2, is provided. This unit contains a tape reader 
for reading process-tape information into the computer and a punch 
for punching out the control tape. The magnetic-tape handler sup- 
plements the internal storage of the computer, greatiy enlarging the 
library of subroutines that can be handled automatically 

The typewriter shown in the foreground is normally furnished 
with the G-15D computer. It contains the operating controls and 
provides a means for typing information out of the computer during 
computation. This latter feature is used, for example, in typing out 
co-ordinates of the cutter as the computer produces the control tape. 
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Figure 6. Photograph of G-15D, AN-2, and MTA-2 


COMPAC-1 PROCESS SHEET 


Figure 7 illustrates the form of the manuscript or process sheet 
that has been developed for use in the Bendix system. The term 
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COMPAC (COMputer Program for Automatic Control) refers gen- | 


erally to the Bendix tape preparation system. COMPAC-1 refers 
specifically to a set of computer routines for the profiling and pock- 
et milling of essentially two dimensional parts whose contours can 
be described in terms of circular arcs and straight line segments. 
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E. C, JOHNSON 


The COMPAC-1 process sheet is divided into five major areas, 


The area in the upper left corner contains basic setup information| 
which is common to the rest of the process data on the sheet. The! 


area below in the left half contains most of the dimensional infor- 


mation as well as certain codes for both profile and pocket milling, | 


The next field contains feed rate and tool diameter; the fourth, 
auxiliary-function commands; and the fifth, information used in 
pocket milling only. 


EIA SAMPLE PART 


Figure 8 illustrates the D-shaped part designed to serve asa 
basis for discussion at the EIA Numerical Control Symposium. As 
indicated by the drawing, a one-inch plate is to be milled down toa 
thickness of three quarters of an inch, leaving only the D-shaped 
area uncut. 
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Drawing of EIA Sample Part 
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Outside Cut 

Consider first the machining of the outside of the D, Figure 9. 
The process planner elected to do this in one pass using a cutter 4.5 
inches in diameter. The point on the x axis at the left numbered 1 
was arbitrarily selected as the setup point, or initial position of the 
tool. The origin of the co-ordinate system might orcinarily have 
been chosen, but in this case the cutter would have interfered with 
the work blank. The co-ordinates of the setup point, x = -1.5 and 
y = 0, are entered in the setup field of the process sheet. The z 
dimension, 0.75, represents the plane of the finished cut, corre- 
sponding to the bottom of the cutter. 
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The column headed clearance plane refers to a plane parallel to 
the x-y plane which is sufficiently far from the workpiece and any 
holding fixtures that motion of the cutter center in this plane is free 
of obstruction. During pocket milling, transitional motions to and 
from the pocket are programmed in this plane automatically by the 
computer. In this example, 1.25 inches provides ample clearance. 
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The next column, tolerance, refers to the mathematical accuracy 
to which curves are to be approximated by straight lines. In this 
case the process planner has called for 0.0002 inch, which will re- 
sult in a very good surface finish. The “1” appearing in the next 
column is a code calling for a chord type of approximation rather 
than tangent or secant. 

The word “print” entered on the first line of the process sheet 
proper indicated to the computer that it is to print out cutter-center 
co-ordinates as it proceeds along the computation. Cutter-center 
points corresponding to all co-ordinates entered on the process 
sheet will then be typed automatically. The “20,” just below, in 
addition calls for every twentieth point along an arc to be printed 
out. The “45” is a code indicating to the computer the significance 
of the “20” — actually the number of the routine to be used in the 
process of computation. 

On the next line over toward the right is entered the feed rate in 
inches per minute and the tool diameter in inches. This information 
is followed by the routine number “40.” The auxiliary function “05” 
called for on the next line causes a reference code to be punched in 
the control tape. In rewinding the tape on the machine control unit, 
it can be made to stop automatically at this point. 

The last five lines shown on the process sheet in Figure 9 
specify points 2 through 6 on the part in the sequence of the cuts to 
be made. (These points have been re-numbered from Figure 8.) 
Point 2 is merely an extension of the base of the D. The computer 
automatically offsets the center of the cutter by the radius so that 
the first motion of the cutter is to the point 2'. Following succes- 
sively around the workpiece, the x, y, and z co-ordinates of points 
3, 4, 5, and 6 are given. The primed numbers represent the cor- 
responding locations of the cutter center, as automatically deter- 
mined by the computer. For the circular section whose end point is 
5, the radius of the circle is aiso entered on the process sheet. The 
minus sign indicates that the circleis to be traversed in the counter- 
clockwise direction. 

In the arc-length column, “0” indicates a straight-line section 
and “1” and arc less than or equal to 180 degrees. An arc greater 
than 180 degrees would be indicated by “2.” The “1” in the next 
column indicates an outside corner. The “2” in the next column 
indicates that the contour of interest lies in the x-y plane, rather 
than the y-z or x-z planes, and that the cutter center is to be offset 
in the x-y plane. The minus sign on the first “2” indicates that this 
is the first point of a continuous sequence. The R number or routine 
number for all information in this field is “00.” This represents 
all of the information required as input to the computer for the ma- 
chining of the outside of the D. 
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Inside Cut 


In machining the inside of the D, use is made of the pocket- 
milling portion of COMPAC-1. As indicated in the first line cf the 
process sheet, Figure 10, the cutter is first changed to one having a 
diameter of one inch in order to generate the small circular seg- 
ments at the base of the D. The minus sign on tool diameter indi- 
cates that the tool is to be offset to the left of the part when looking 
along the direction of the cut. In the previous operation of ma- 
chining the outside of the D, the cutter must be offset to the right. 
For that case then the tool diameter carries a positive sign, as 
indicated in Figure 9. 
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Figure 10, EIA Sample Part, Inside Pocket 

The preceding sequence of cuts left the center of the cutter at 
point 6'. A machine stop indicated by “04” in the auxiliary function 
column is now called for in order to make the necessary cutter 
change. The pocket-milling sequence proper is begun by entering a 
line of information in the extreme right field of the process sheet. 
The depth of cut on the final pass around the pocket was arbitrarily 
chosen to be 0.25 inch. Because of the quarter-inch bottom radius 
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of the cutter, at least three roughing passes are then required to 
completely clean out the bottom of the pocket. The number of sec-| 
tions making up the boundary of the pocket is four: one straight | 
line and three circular sections. | 

These four sections are defined in the next four lines of the 
process sheet just as they would be for profile milling. The x, y, 
and z co-ordinates of the points 7, 8, 9, and 10 on the desired fin- 
ished surface are listed in that sequence. Circle radius and the 
coded information are treated just as they would be in profile mil- | 
ling. In addition, the initial-clearance column is filled in for each 
section in the pocket boundary. The initial clearance, as indicated 
in the diagram, is in effect the distance between the part surface and 
the center of the cutter on the first or inside pass. The exact value 
given for initial clearance is not critical, the essential requirement 
being that the first rough cut clean out the center of the pocket. 

In this example the entire pocket can be completely machined | 
ignoring the small circular sections since they are actually formed 
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by the shape of the cutter itself. An initial clearance equal to the | 
radius of the cutter on these two sections causes them to be skipped | 


automatically by the computer. 

Given this amount of information — consisting essentially of the 
description of the pocket boundary, the number of roughing cuts to 
be made, and the depth of the finishing cut —the computer automati- 
cally generates all of the motions required to machine the pocket. 
Starting from point 6', the cutter is first withdrawn to the clearance 
plane. It then moves directly to a point near the center of the pock- 
et. The plunge into the pocket to the required depth takes place ata 
feed rate which is automatically reduced to 25 percent of the value 
entered on the process sheet. The cutter then follows the sequence 
of motions indicated to complete the three roughing cuts. In ad- 
vancing to the final cut, the cutter slows down again to 25 percent of 
the specified feed rate. It then proceeds around the required path. 

In order to avoid an overshoot which might damage the work- 
piece at the two corners, automatic slowdowns are programmed by 
placing minus signs before the corner codes on points 7 and 9. The 
computer then automatically slows the cutter down to a feed rate of 
2 inches per minute, beginning at a distance of one eighth of an inch 
before the corner. These slowdowns are ignored on all the roughing 
passes, being taken into account only on the final cut. Likewise, the 
approximation tolerance on the arcs is automatically increased to 
one eighth of the cutter radius on all except the last roughing cut 
and the final cut, where the specified tolerance is used. The number 
of points that have to be determined on circles involved in the early 
roughing cuts is thereby reduced. 

The pocket-milling sequence ends with an automatic withdrawal 


— 





to th 
meré 
com 
the v 
mati 
tatio 
V 
gene 
The 
boun 
addit 
natic 
is ar 
] 
for < 
The 
the 
tape 
wert 
mor 
sam 


by t 
punc¢ 
with 
vide 
cher 


plan 
15 | 
min 
pute 
of t 
the 
cen 
part 
pro 
41 f 


it is 


ment 


1ined 
rmed 
o the 
pped 


f the 
its to 
nati- 
cket. 
‘ance 
OCK - 
pata 
value 
lence 
| ad- 
nt of 
ith. 
rork- 
d by 
The 
ite of 
inch 


ghing | 


», the 
ed to 
z cut 
m ber 
early 


‘awal 


DATA PROCESSING FOR NUMERICAL MACHINING 15 


to the clearance plane. The last two lines on the process sheet 
merely cail for motion back to the original setup point so as to 
complete the circuit. The process sheet is terminated by entering 
the word “end” at the bottom. This code causes the computer auto- 
matically to cease reading the process tape, complete the compu- 
tations, and punch a machine stop code on the control tape. 

With just this amount of information, the computer is able to 
generate the motions required to completely machine the pocket. 
The only precise information needed is that concerning points on the 
boundary of the pocket which were given in the original drawing. No 
additional manual computation is necessary except for the determi- 
nation of suitable initial clearances. As mentioned previously, this 
is an operation in which no precision is necessary. 

Figure 11 shows the appearance of the finished process sheet 
for all of the operations involved in machining the D-shaped part. 
The information required fits on one sheet. The typewritten copy of 
the process sheet which results during the punching of the process 
tape on the Flexowriter is shown in Figure 12. The column headings 
were added in a subsequent operation merely to make the sheet 
more nearly self-explanatory. The numerical data is exactly the 
same as that appearing on the handwritten process sheet. 

Figure 13 shows the record of cutter-center locations produced 
by the computer as it goes through the process of computing and 
punching the control tape. This record was produced in accordance 
with the “print” instruction entered on the process sheet. It pro- 
vides a means for monitoring the operation of the computer and 
checking roughly on the accuracy of the finished control tape. 

The total amount of time required on the part of the process 
planner to complete the handwritten process sheet for this part was 
15 minutes. The time required to punch the process tape was 5 
minutes, and the time required to verify it, 7 minutes. The com- 
puter required a total time of 37 minutes, which includes the read-in 
of the process tape, performance of all the calculations involved, 
the punching of the control tape, and the print-out of the cutter- 
center locations. Total time to produce the control tape from the 
part drawing as a start was therefore 64 minutes. The length of the 
process tape was 6 feet, and the total length of the control tape was 
41 feet. These statistics are summarized in Figure 14. 


SUMMARY 
The Bendix tape preparation system is extremely flexible since 


it is based on the use of a stored-program, general-purpose, digital 
computer. Such a computer can be programmed to handle virtually 
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any computation of practical significance. The term COMPAC-1, as 

already explained, applies toa particular set of computer routines 

developed for essentially two-dimensional part contours that can be 

described in terms of straight lines and circular arcs. With ad- 

ditional programming, the system can be readily extended to include 

more complex two-dimensional parts, as well as three-dimensiona! | 
work, 

Experience to date with COMPAC has indicated that the system 
does represent an effective method of producing tapes for numerical 
control. It is expected to become still more effective as additional | 
programming is completed. 
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SUMMARY OF TAPE PREPARATION 




















PROCESS PLANNING AND PREPARATION 
OF HANDWRITTEN PROCESS SHEET 15 MIN. 
PUNCH PROCESS TAPE 5 MIN 
VERIFY PROCESS TAPE 7 MIN 
COMPUTE AND PUNCH CONTROL TAPE 37 MIN 
TOTAL TIME 64 MIN 
LENGTH OF PROCESS TAPE 6 FT 
LENGTH OF CONTROL TAPE 41 FT 
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Cam Design on a Digital Computer 


By Walter Hoffman 
Department of Mathematics 
Wayne State University 


The determination of cam outlines generally involves long and 
tedious calculations which are particularly well suited to calculation 
on an electronic digital computer. 

The particular application to be reported in this paper is con- 
cerned with the production of a cam with previously determined lift 
characteristics on acam grinding attachment which is shown in Fig- 
ure 1 in somewhat simplified form. The master cam and its shaft 
are axially connected to the cam shaft which is to be ground. The 
grinding wheel is positioned against the cam shaft and the grinding 
wheel motor is started. The cam shaft drive is then started and 
operation begins. Since the roller follower rotates on a fixed axis, 
the master cam and the cam shaft pivot away from their original 
position when the cam contacts the follower on the rise of the cam. 
A strong spring force keeps the master cam in contact with the 
roller follower, and thus controls the contour which is ground on the 
design cam. 

The problem to be solved is the determination of the master cam 
outline for a given design cam contour. The obvious complications 
are the difference in size of the roller follower and design cam, as 
well as the lack of alignment of their axes. An additional compli- 
cation is the motion of the pivot arm, which produces a phase shift 
between master cam and design cam. When the nose of the design 
cam is in contact with the grinding wheel, the master cam has 
turned approximately 1 1/2 degrees past the nose against the roller 
follower. This phase shift changes continuously, thus complicating 
the computations, and obviating the possibility of obtaining even 
moderate accuracy by means of conventional layout procedures. 

Figure 2 represents the design cam in contact with the grinding 
wheel. We are given 


Ro = radius of grinding wheel, 


L, (9) = design cam lift as a function of angular displacement 
at one degree intervals, 


180 d 
e (0) - or - - cam velocity, given at one degree intervals. 


For example, in Figure 2, e (@) = RS 
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STMALIFILD SKETCH OF GRUMLING ATTACHMENT 
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Figure 1. Simplified Sketch of Grinding Attachment. 
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CAM 


We first compute 


e (0) 
y, (0) = arctan ‘ 
| RG + Lo(@) 
which enables us to compute the center to center distances 


Br, (8) = (Ro + (8)) sec 9 


Since 9-9, =¥W, 2,,, is now tabulated as a function of ¥. 
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ign Cam in Contact with Grinding Wheel. 
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A DIGITAL COMPUTER 3 

Figure 3 depicts the front view of the mechanism shown in Fig- 
ure 1. The center distance 2,,, has been subjected to an angular 
displacement , which is due to the heretofore neglected pivotal 
motion of the camshaft. A rectangular Cartesian coordinate system 
with origin at the center of the grinding wheel is used to compute 
the path of the moving camshaft axis. Its coordinates x (y) and y (W) 
are seen to be a solution of the system 


where 


é Clearly, the desired solution of this system is that for which y 
: is non-negative. 
We now compute a (yw), the center to center distance of master 
camshaft and roller follower as 
2 (Wy) jy” + (x - c)* 
.— where 


S Ro + Re = Ry - Ry 
‘ The angular displacement py, of a is next computed as 


‘ y 
b, ei arctan — 
x 


and 2 (W) is tabulated as a function of YW + u, = P. 
Figure 4 illustrates the master cam in contact with its roller 
follower. Essentially the initial steps have to be reversed, which 


; oa a 
requires the computation of the master cam velocity ,,. This is 


Z 
dp° 
accomplished by Lagrangian interpolation. The remaining steps are 
to calculate the quantities 


L; (P) =a cos¥,- Rr. 


Finally, L; is tabulated as a function of P+, , and a last inter- 
polation is performed to compute a list of values of L; at successive 
integer values of P + ¥, , which are the required results. 

On the basis of a few hand calculations for checking purposes, 
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Figure 4, Master Cam in Contact with Roller F 


TER 


the time required to solve this problem on a desk calculator was 
estimated in excess of seven months. The problem was coded for 
the IBM 650 EDPM in fixed point arithmetic and required approxi- 
mately thirteen minutes of running time. 
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What the Test Engineer Needs in the Way of Statistics 


By Leonard G. Johnson 
Research Staff 
General Motors Corporation, Detroit 


INTRODUCTION 


Experimental science is concerned with drawing generai conclu- 
sions from particular observations. This process of going from the 
particular to the general is what is known as induction. Hence, any 
science in which general statements are formulated by inference 
from observations is an inductive science. 

The engineer who subjects specimens to an endurance test in 
order to determine their fatigue life is engaged in an inductive 
science. He makes particular observations as to the hours or cycles 
to failure of particular specimens under a given set of conditions, 
and from this data he draws general conclusions as to the distri- 
bution of life, or more specifically he estimates such quantities as 
the average life, or the median life, of the entire population from 
which these sample observations supposedly originate. 

The process of arriving at population estimates from the data 
furnished by means of samples falls within the branch of applied 
mathematics known as Statistics. Hence, statistics makes use of 
inductive reasoning. The planning and interpretation of fatigue tests 
is an engineering application of statistical inference, and conse- 
quently the test engineer needs the help of certain concepts in 
mathematical statistics. What some of the necessities are is the 
subject of this report. 


STATISTICAL ESTIMATION VIA SAMPLES 


The first statistical need of the test engineer is a method of 
estimating the population of a set of fatigue life observations ob- 
tained under fixed conditions of material, loading, speed, tempera- 
ture, humidity, lubrication, etc. Strictly speaking, determining the 
exact nature of a fatigue life population defies all human ingenuity. 
All that we are ever able to do is to arrive at a convenient mathe- 
matical model based on certain assumptions concerning fatigue 
phenomena. This mathematical model will contain certain constants 
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LEONARD G, JOHNSO} WHAT 
or parameters which are estimated from the set of observed fatigu® patur: 
lives. Hence, the real problem of population estimation is glosseiR what 
over by means of assumptions which lead to special types of equa- to be 
tions for fatigue life distributions and then the estimation problem T 
is reduced to merely an estimation of parameters in a mathematical lated 
formula. The validity of the final estimate is therefore contingent 


hit 
upon the validity of the original assumptions concerning the nature . | 
of the population. One of the simplest assumptions which can be the s 
made about fatigue phenomena is that a specimen consists of af assur 
series of elementary particles or “links” held together to formaf} jives 
“chain.” Since a chain is no stronger than its weakest link, it is longe 
reasonable to assume that every fatigue failure is initiated by the of po 
failure of the weakest link. Hence, it is assumed that fatigue failuresB tract 


are examples of extreme values, i.e., they are smallest strength of the 
values. This assumption leads to the Weibull distribution as an a the 
asymptotic limit, reference 1. According to this particular theory, if 





sci 

we assume that the minimum possibie life isa, then the probability a 
of failure within time x >a (which may be measured in hours, stress too }, 
cycles, or some other unit proportional to time) is given by F(x,) 
q > value 

F(x) 1 - Exp | - ( r) ) } (1) any 

posit 

In general, all three of the parameters a, b, and 0 of the Weibull V 
equation (1) are unknown and the test engineer must resort to actual fatig 
test results to estimate them. In the case of specimens for whicha grap 


very short life, arbitrarily close to zero, is observed, we may take ] 
a = U0, in which case (1) becomes 


men. 
F(x) = 1 - Exp[-(x/@)”] (x>0) (2) 
What the test engineer likes to have is a graphical method of esti- 
mating the parameters b and @ in this formula. (2) can be put into ; 
the form 
log. log, 1-Fi&) c + b log, x, (3) 
The 
where c = - b log,9. five 


Thus, a Weibull distribution (2) is equivalent to a linear relation 
between the logarithm of fatigue life x and the logarithm of the 
1 
logarith ———— aa - Te sey sar a 
—_ = « 1-fraction failed at x . en eo 
constructed on this basis, i.e., by making the abscissa a logarithmic 
scale representing log.x, and by making the ordinate a log log scale 


1 ; : 
representing log, lOBe TTF(x) * At this point the question which 
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naturally arises is, “If a fatigue life x is observed, how do we know 
what value to assign to F(x) in (3), since F(x) is the unknown function 
to be determined?” 

The answer to this question is that if we are given only one iso- 
lated value of fatigue life x out of a sample of N values, we can say 
nothing about the fraction of the population failed within time x un- 
less we know the relative order or rank of the fatigue life x within 
the sample of N observed values of which it is a member. Let us 
assume that for a random sample of size N, the observed fatigue 
lives in order of length of life, i.e., from the shortest life to the 
longest life, are x,, X,, X,...., Xn- It can be shown that regardless 
of population, there is a unique median value for F(x ), i.e., for the 
fraction of the population having a life less than or equal to the life 
of the first failure in N. This means that when we assign to x, such 
atheoretical median value of F(x,) we are just as likely to give the 
abscissa x, too high an ordinate as too low an ordinate on Weibull 
probability paper. Hence, with 1: 1 odds of being too high vs. being 
too low in the population, we accept this theoretical median value of 
F(x,) as a reasonable plotting position for fatigue life x,. Median 
values have also been determined for F(x,), F(x,), ...., F(xn) for 
any sample of finite size N, regardless of population. These plotting 
positions are derived in reference 2. 

With the median plotting positions known for a set of observed 
fatigue lives x,, X,, X;,..-.,%n, it is a simple matter to make 
graphical estimates of the parameters b and 0 of (2), 

For example, let us assume that a random sample of five speci- 
mens yields the following observed fatigue lives: 


X, 51 hours 
x, = 97 hours 
x, = 150 hours Abscissas on Weibull Probability Paper 
x, = 220 hours 
x, = 300 hours 


The median plotting positions (median ranks) to be assigned to these 
five values are: 


F(x,) = .1294 
F(x,) = .3147 
F(x,) = .5000 Ordinates on Weibull Probability Paper 
F(x,) = .6853 


F(x,) = .8706 (See Table 1) 











TABLE 1 — MEDIAN RANKS 


Sample Size =n 
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onour wn = 


19 


1 2 3 4 5 6 7 8 9 10 
1 .5000 .2929 .2063 .1591 .1294 .1091 .0943 .0830 .0741 .0670 
2 7071 .5000 .3864 .3147 .2655 .2295 .2021 .1806 .1639 
3 .7937 .6136 .5000 .4218 .3648 .3213 .2871 .2594 
4 .8409 .6853 .5782 .5000 .4404 .3935 .3557 
5 8706 .7345 .6352 .5596 .5000 .4519 
J 6 8909 .7705 .6787 .6065 .5481 
7 9057 .7979 .7129 .6443 
8 9170 .8194 .7406 
9 9259 .8368 
10 9330 
Sample Size = n 
oe | 6 oe 
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The actual plot of these abscissas and ordinates on Weibull 
probability paper appears in Figure 1. The straight line represents 
an estimate of the population from which the data supposedly origi- 
nates. By reading the graphic picture in the Figure, we can predict 
what percentage of bearings will fail within any specified number of 
hours. For example, we predict 32 per cent will have failed within 


100 hours. 
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We note that the plot exhibits approximate linearity, which fact 
bolsters our confidence in the assumption of a Weibull distribution. 
The slope of the line is an estimate of the parameter b in(2). Hence, 
the parameter b has been given the name Weibull slope. In order to 
estimate the parameter 9 in (2), we note that if we put x =9@ in (2) 
we obtain F(Q) = 1 - e~* = 0.632. Thus, @ is the abscissa of the 
the Weibull plot which has an ordinate equal to 63.2 per cent failed. 
It has come to be known as the characteristic life. 

Furthermore, the median life is the abscissa on the Weibull line 
which corresponds to the ordinate 50 per cent failed. This graphical 
method of estimating the parameters of a fatigue life distribution 
enables the engineer to get a feeling for the statistical concepts in- 
volved in this process of induction. No matter what method is used 
to fit a straight line to the plotted points, the engineer sees the un- 
certainties inherent in statistical estimation, especially if he carries 
out the procedure for numerous samples from the same population. 
The observed slope and median life are never constant, even for 
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samples taken from the same population. This raises the next 
question for which the test engineer needs an answer, namely, “How 
close to the true population values are the particular parameter 
values which we happen to observe for a sample of size N?” Ip 
other words, how wide a confidence band about an observed Weibull 
line is required in order that we can have 90 per cent probability, or 
any other specified probability, that the true population abscissa cor- 
responding to any ordinate will fall within the band so constructed? 
If it is assumed that the Weibull slope is known exactly, this question 
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is answered by the type of charts shown in Figures 2 through 4 for 
Weibull slopes 1.0, 1.5, and 2.0, respectively. In the general case, 
where the Weibull slope is estimated by the sample at hand, and 
hence subject to error, the confidence band will be still wider be- 
cause the population could have a smaller Weibull slope than that 
observed. The test engineer needs to know what this total uncer- 
tainty range is for a sample of any finite size N. Once these un- 
certainty ranges have been determined for all values of N, the 
engineer is able to determine what size sample is required for any 
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fatigue test in order that the total error in the estimation of am 
population life abscissa at a given ordinate (per cent failed) does ng 
exceed some tolerable value. 


SIGNIFICANT DIFFERENCES 


A probiem of prime importance in the field of fatigue life testing 
is the development of materials of improved endurance. This brings 
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of any® another statistical question into focus, namely, that of determining 
es no how confidently we can state that an observed improvement indicates 
a veal improvement. It must be admitted that two samples from the 

same population can show quite a difference in mean life. But this 

» difference is only due to random variations within one population, 

| and does not arise due to population differences. In searching for 

) betier materials we are really looking for better populations. 
esting Hence, whenever the test engineer observes that a new material 
br ings seems to be superior to some standard material, he needs to know 
how great the observed improvement must be in order that the 
probability of its occurrence by luck, i.e., by pure chance from a 
) population of standard quality, is so remote that it is more reason- 
able to assume that a better population, i.e., a better material, has 
been found. This question is partially answered by the type of chart 
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shown in Figures 5 through 10, in which exact Weibull slopes are 
assumed to be known in the range 1.0 to 2.0. All these charts as- 
sume that both samples under comparison have the same known 
Weibull slope in the interval 1.0 to 2.0, as well as a zero minimum 
life. The source of these charts is reference 3. These charts en- 
able the engineer to estimate how many specimens he must have in 
the samples being compared in order to detect an improvement of 
specified magnitude with 90 per cent confidence, or whatever hap- 
pens to be the acceptable degree of confidence. In the most general 
case only estimates of the true Weibull slopes are obtainable, and 
since these are subject to error, extra allowance should be made by 
using a chart corresponding to a smaller Weibull slope than that 
observed for either sample under comparison. 

Having computed the mean life ratio, we can refer to Figures 8 
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through 10 to obtain the confidence number corresponding to this 
mean life ratio. (This confidence number is the probability that the 
10wn | true mean life ratio of the populations is greater than unity.) 


num 

en- ’ 

re in CONC LUSION 
it of 
hap- We have indicated some oi the statistical problems involved in 
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For engineering purposes we have proposed the following methods 
in dealing with these problems: 


1. Estimate parameters graphically on Weibull probability 
paper. In case the minimum life is greater than zero, try sub- 
tracting off various amounts from the observed fatigue lives in 4 
sample until a linear Weibull plot is obtained. The subtrahend 
yielding the most nearly linear tendency is then an estimate of the 
minimum life a in equation (1). 


2. Determine upper and lower limits on any quantile level by 
employing confidence bands of the type shown in Figures 2 through 
4 with proper allowance for errors in Weibull slope. 


3. Determine the confidence assignable to an observed im- 
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Figure 8. Weibull Slope 1.6 
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provement by using charts of the type shown in Figures 5 through 
10, employing a lower limit on Weibull slope. 


It should be borne in mind that these are only suggestions which 
serve aS a guide for the test engineer. A rigorous and compre- 
hensive mathematical treatment of these problems is beyond the 
realm of this report. 
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Error Correction in Solutions of Linear Systems 


By James S. Seward 
Engineering Staff 
Ford Motor Company 


Introduction 


After solving a large-scale system of linear equations on a 
computer, it is frequently discovered that an error was made in 
calculating or recording a coefficient of one of the unknowns. A 
single error usually affects the entire solution and the problem has 
to be recomputed to obtain correct answers. A considerable amount 
of machine time and money can be saved by using techniques which 
correct a solution to a system of equations directly without re- 
solving the system. The techniques also have application to studies 
todetermine the influence of a particular coefficient or coefficients 
on the solution of a linear system. A system can be solved using an 
approximate value of a coefficient. An error-correcting technique 
can then be used to evaluate the effects of a change in a coefficient 
without having to re-solve the system by direct methods. 

Two error-correcting techniques are discussed here. The first 
method, which is easily carried out using a desk calculator, corrects 
the solution of a system of linear simultaneous equations to corre- 
spond to a single change made to a coefficient. The second method, 
best suited for multiple corrections using a digital computer, cor- 
rects the inverse of the matrix of coefficients. The corrected in- 
verse is then used to obtain a new solution by multiplying the new 
inverse by the column matrix of constants. 

The second method has been programmed for an IBM 650 com- 
puter equipped with an automatic floating decimal device and in- 
dexing accumulators. A brief description of the program is given. 


L Correction of Single Errors 
Suppose the system of linear equations that has been solved is 
represented in the usual fashion as: 
AX =K (1) 


where A is an n by n matrix whose elements are the coefficients of 
then unknowns, x,; X is a column vector, x, , X2.... X,; and K is 
the column vector composed of the constants k,, k2... Ky. 
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The system that should have been solved is 
(A +B)X =K (2) 


where B is a matrix composed of all zeros with the exception of the 
one coefficient to be corrected. 

The solution we have obtained to (1) is the column vector 
X ,= A™' K whereas the solution we wish to obtain is 


X =(A+B)" K (3 
Following Dwyer (1), (A + B)~ can be expanded as 
(A + B)™ = (A+ AA™ B)” 
(A(I + A™ B))” (4 
(I+ A™ B)” A” 


where I is the identity matrix. 

For simplicity in the proof, assume that the element to be cor- 
rected is located in position n, n, of the matrix A. The correction 
matrix B then consists of all zeros with the exception of the number 
in the last row and the last column. 

The multiplication of A~ by B will then produce a matrix con- 
sisting of all zero elements with the exception of the last column 
which will consist of the values a>, b,,,,.-- am Dyn where the a; 
are elements of A™ and the b,,,, is the correction value. 

When the identity matrix is added 
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If the values ain bay, , Aon Dnn--- ary n Dnn are combined intoa 
column vector D,-the expression can be written as 
, 
(I+ A™ B) (5 
0 1 + — Baw 


The inverse of this matrix can be found in the following way: 
Let 
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where R, S, T and U are of the same dimensions asI,_,, D, , and 


the sinvle constant (1 - a>. b,.). Then since AA~ must equal I, 
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where s = 1+ a7, Dann and the zero matrices in the identity are of 
the proper dimensions. Then 


LR oo DP @f,., 


I,-, 8 + DU =0 
sT 0 
sU l 
Therefore 
T=-0, 1 ae. Bek S=- : D 
and 
/ 1 
- p \ 
(I- A” Bj)” 1 ) (6) 
0 ; 
Ss / 


1 
The elements of - - D are now of the form 
5s 
where i =1,2,..., n-l 


Recapitulating 
X = (A + B)™ K = (I+ A™ B)™ (A™ K} = (1+ A™ B)” X, (7) 


If the column vector of the n initial solutions, X,, is split into 
wo vectors, one with all the elements except xX, ,, the other with 
fon Only, the system can be written 


[ Xp-, \ [Nn <“ 
\xn / \o 


Xn 
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and multiplied out to give: 
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The need for handling the nth solution x, , in a different manner 
than the other, can be removed by noting that 


Xo on Xo .n -1 
Xn . Ss Xo - “ (a nn Dnn) 


The new values x, can then be given in terms of the old value x 
as follows: 


0.7 


Xo n 
X + Xo.r - = 





qn ’ 
fa. Onn)  =1,..- 0) 

By examining the proof, it can be seen that the restriction on the 
position of the element to be corrected was done for convenience 
and is not a necessary condition. By rearranging the submatrices, 
tne proof can be repeated for an element located in any position. 

The general rule for correcting a single element located in 
position jk of the original matrix of coefficients given (1) the ele- 
ments of the inverse matrix, (2) the original solution vector X 9, and 
(3) the value of the correction b;, , can be stated as follows: 

(1) Compute the correcting factor 

(Xo) (b ji) (Xo) (Dik) 


- 1+art bi, 
(2) The new values of the solution are given by: 


eo Oo. 


» gl 
KX, = Xgr- Car; 


Note that the elements of the inverse that are needed are those 
in the column corresponding to the row in which the correction was 
made. 

The only restriction on the correction factor, b;, , is that it can- 
not equal -1/ajz . If this would happen, s would equal zero and 
|1+ A Bl =0 since one row would be all zero. Then [{A| |I + A™ B 

0 or |A + B| = 0 which indicates that the corrected system could 
have no solution. 


Numerical Examples 
Example 1: 


Given the system: 


2 -10 15 32\ /x,\ / 23 \ 
19 45 -14 -8 X, 57 
-12 16 27 13 XK; 47 
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x, = 2.0 
x,= 1.0 
x, = 3.0 
xX, = -2.0 


and the inverse matrix has the elements 


/ 0.02873 
| -0.00695 
\ 0.01825 
\ -0.00282 


\ 


Suppose the element a,, 


s=l+atib 


c 


0.02436 
0.01239 
0.01440 
-0.02267 


-0.02302  -0.01519 
0.01572 0.00419 
0.00791 -0.02041 
0.01991 0.02322 / 


which is now 16 should have had the value of 
10. The correction to be applied is the addition of -6, i.e., b,, = -6. 
The correction is made as follows: 


k 


=. 2 
0,K 


1 + (0.01572) (-6.0) = .90568 


Dix = -6.62486 


The new values of x are then given by 


X4 


Xo. 
Xo 2 
Xo\3 


Xo 4 


cals 1.8475 
ca, = 1.1041 
ca;, = 3.0524 
ca,, = -1.8681 


These results agree to 4 decimal places with the solution of the 
modified system obtained by conventional means. 


Example 2: Change of Sign 


A frequent error in the preparation of data is recording the 
wrong sign of the number. 
minus twice the value of the original coefficient. Given the system 
inexample 1, suppose the element a,, which was given as +16 should 


have been -16. For this case: 


Ds. 


c 


-32 s 
Xo.2 
—a Dos 


In this case the correction factor is 


1 - 32(a35) = 0.49696 


32 
649696 -64.39150 
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Then 

K,=Xo.- aj, = 0.5177 
Xo = Xo,2- Caz, 2.0122 
X3=Xo,3- Cayz = 3.5093 
Kq=Xo4- Cas = -0.7180 


These results can be verified by substitution into the modified 
system of equations. 


Il. Correction of Multiple Errors 





The method presented for correction of a solution of a system of 
linear equations cannot be used more than once on a system since it 
does not correct all the elements of the inverse matrix to corre- 
spond to the new coefficient matrix. To correct for multiple errors, 
the inverse matrix must be modified and a matrix multiplication 
carried out to give the corrected solution. 

From the previous development for a correction in position nn, 


(A +B)" =(I1+A™" B)“* A‘ 


= - - p| A” 
Ss 
1 
\9 = 


Partitioning A™ intofour submatrices E, , E, , E, and E, , where 
E, is size n-1 by n-1, E, is n-1 by 1, E, is 1 by n-l and E, az,,, the 
product above can be written 
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; DE , ‘ 
the correction - = 3 which is to be subtracted from the inverse 


submatrix E, is 





qu! a=! 5 71 =! cl i=l 
ain 4ni Onn Ain Anz Pnn +++ Ain @n n-: Don 
1 on 1 ay, t i +a, 2 
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The correction for the inverse submatrix E, has elements of the 
same form 





“1 aol » 
ain ann ann 
-1 
1 + ann Dan 
. DE, J 
: 
\ an at 


\ . . 
\ an-in 4nn Din 





-1 
\ 1 + ann Dan 


The correction for E, and E, consists in merely dividing the 
elements of this row of the inverse by s =1 + —. b,,- So that a 
different type of correction does not have to be applied to this row 
the elements of E, which are of the form, 

an k 
—_———— (k=1,...n-1), can be rewritten as 
1+ ann Dan 


=1l -l " 
ann 4nk Pnn 


=} 


- all 
“_ 1+ ann Pnn 


a 


Thus, for a correction b added toa,, , the element in the n,n po- 
sition of a matrix A, every element, a, ; of the original inverse must 
have the factor 


pe 
71 ol t 

. 4pn 4ng Pnn 
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added to it. 
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If the element that is corrected is in the j,k position of the origi- 
nal matrix of coefficients, the correction factor is 


baal 5 


aw 
7 an} Akg Dix 


= (p,q = 1... a) 
1 + ay | Diy 


This is the formula developed by Sherman and Morrison in (2) 
using summation methods. Since the inverse is corrected each time 
a correction is made to the matrix elements, the technique may be 
used any number of times on the same system. 


Numerical Example 


Given the matrix of example 1 and its inverse, suppose both the 
elements in position 1,2 and 3,4 are to be corrected. Let b,, = 10 
and bs, = -3. The coefficient matrix after both corrections is then 


26 0 15 32 
19 45 -14 -8 
-12 16 27 10 
32 29 -35 28 


To correct for the change in position 1,2, the correction factors 





are 
© « ¢ oe 
a ap ,1 42,q Dis (p,q 1,2 eee n) 
Zé ar Dio 
- (10.7469)a5\ a3; 
Then c,, -(10.7649)a;\, a3; = 0.00215 and the new value for aj, is 


0.02873 + 0.00215 - 0.03088. The other elements are corrected in 
the same fashion to give: 


0.03088 0.02053 -0.02787 -0.01648 
-0.00747 0.01331 0.01689 0.00450 
0.01961 0.01197 0.00482 -0.02123 
-0.00303 -0.02229 0.02039 -0.02334 


Al 


The correction for the second change b,, = -3 can now be ap- 
o ’ 
plied using this newly-computed inverse. 
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al o=l 1 
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1 +a b,, 
= (3.19547)a5) az, 
Thus c,, = (3.10547)aj, a3, = 0.00027 and aj, becomes 0.03088 + 
00027 = 0.03115. 
The inverse after both corrections is 

/ 0.03115 0.02252 -0.02969 -0.01856 

-0.00764 0.01211 0.01799 0.00576 

s 0.01957 0.01163 0.00514  -0.02087 

\ -0.00323 0.02375 0.02171 0.02487 


Since the method of correcting inverses is an exact one, the 
results agree with the results of inverting the corrected matrix to 
within the round-off error of each method. 

The new solution of the system of linear equations can now be 
found by carrying out the multiplication of the constant vector, K, by 
the corrected inverse. The final solution is 


X) 1.8667 
X2= 0.9684 
Xs 2.7735 
X4 = -2.0980 


Extension of the Method 


One interesting aspect of this method for the correction of 
multiple errors is that it suggests a different approach to the in- 
version of matrices. Given a square matrix of order n, an identity 
matrix could be assumed to be the initial inverse. The inversion 
process could be accomplished by reasoning that the original matrix 
was not the identity and applying the n* corrections to adjust the 
inverse for the actual value of the matrix. Using the correction 
technique as outlined here on the order of n* multiplications would 
be required compared to n° for the conventional Gauss-Jordan 
method. However there are some advantages for the inversion of 
Sparse matrices or the inversion of systems similar to some whose 
inverses have already been obtained. The method does not require 
the use of the identity for a starting point. Any known matrix- 
inverse pair could be used. 
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The Computer Program 


The present program for solving systems of linear equations on 
an IBM 650 computer with floating decimal arithmetic and indexing 
accumulators has been designed so that the inverse matrix which is 
developed in the process of solving the system can be punched out 
on cards if the console switch is set to minus. These output cards 
are of the same form as the input cards to the linear equation 
solver, specifying the row and column number of an element and its 
value in floating-decimal form. It is a good practice to specify that 
the inverse matrix be punched and saved for every large linear 
system. 

To operate on an inverse for either, (1) correcting for errors in 
the original inverse, or (2) studying the effect of varying values in 
the matrix of coefficients, a card is prepared for each change to be 
made in the original system specifying the position of the element to 
be corrected and its value in floating-decimal form. [If a new so- 
lution corresponding to the inverse matrix is desired, an indicator 
is punched on the input card. The original inverse and the constant 
vector, K, are read into the computer. The program then reads the 
first correction card and corrects all of the elements of the inverse 
accordingly. If a new solution is desired, the machine then multi- 
plies the new inverse by the constant vector and punches out the 
solution. The next correction is read in and carried out in the same 
manner. The last correction card should always have the indicator 
punched so that the final solution vector will be punched out. 


Time Estimates 


The time required to solve atest system of 38 equations in 38 
unknowns is approximately 35 minutes on the 650 using the present 
program which works from the largest pivot element down. To cor- 
rect one error ina 38 by 38 inverse using the correction program 
requires about 55 seconds. Thus as long as fewer than 38 elements 
are to be corrected in the system, it is more economical of machine 
time to use the error-correcting program than to re-solve com- 
pletely the corrected linear system. It can be seen that the number 
of corrections that can be made economically with the present pro- 
gram is equal to n, the order of the system. 
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The Automobile Engine as a Dynamic Vibration 
Absorber ' 


By John A. Carlson 

Chief of Engineering Analysis Department 
Research and Development 

Teletype Corporation, Chicago 


ABSTRACT 


In all fields of engineering modern computational techniques 
have extended the usefulness of mathematical analysis in design. 
The following paper presents an example of how an analytical study 
can solve a design problem and clarify the behavior of a complicated 
mechanical system. It is shown that the automobile engine mass on 
its flexible mountings can be used as a dynamic vibration absorber 
to reduce shake of the vehicle frame caused by a mode of vibration 
of the car wheels on the tires known as wheel bounce. The conclu- 
sion is reached that the same analytical techniques could be used to 
study the absorption of other modes of vibration by the engine, as 
well as many other vibration problems in vehicles. 


INTRODUCTION 


In a previous paper by the author (1) it was found that motions of 
the vehicle could be obtained with satisfactory accuracy for design 
by representing the vehicle by relatively simple mechanical systems 
and applying forcing functions to these systems. Low frequency 
motions of the frame an body were reproduced and forces on the 
frame caused by motions of the unsprung masses at higher frequen- 
cies were recorded through analysis using a computer. The vibra- 
tions of the unsprung masses were also found to be accurately rep- 
resented. Changing certain parameters in these simple mechanical 
systems was found to have significant effects on the resultant 
motions. The conclusion was reached that the results of analyzing 


This paper is based onresults presented ina thesis submitted in 
May, 1955, in partial fulfillment of the requirements for the degree 
of Doctor of Philosophy at the California Institute of Technology. 

(1) Numbers in parenthesis refer to similarly numbered refer- 


ences in bibliography at end of paper. 
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the vehicle were of sufficient accuracy that design studies could be 
performed on computers and the findings applied directly to the 
vehicle. An example of such a design study is presented here, 
demonstrating the application of analysis to a specific problem. 

For low frequency body motions and motions of the unsprung 
masses in response to road excitation, the assumption that the 
frame, body, and other components mounted on the frame and in the 
car constituted a rigid mass was found to hold. Of course, actually 
it may be expected that other vibrations exist in these components 
and that some of these vibrations may influence the motions of the 
frame. 

The engine on its mountings is a mass-spring system of sizable 
proportions and its motion might be expected to have a noticeable 
effect on the motion of the frame. Some experimental work has been 
reported (2) indicating that this mass-spring system might be used as 
a dynamic vibration absorber. Recently, some work has been done 
on tuning engine mounts (3) in which some of the sources of vibration 
are cited. In the present study, small! vertical motions of the engine 
on its mounts, arising from the vertical movement of the chassis 
caused by road excitation, was considered. The vibrations caused 
by the running of the engine were neglected. Symmetric motion of 
the vehicle was studied. 


ANALYSIS OF THE VEHICLE 


Figure 1 shows a schematic diagram of the automobile as a five- 
degree-of-freedom mechanical system. Here, the engine mounts 
were replaced by a single spring and a damper, K, representing the 
combined vertical spring rate and G, the combined vertical damping 
of the mounts acting in parallel; M, is the mass of the engine and 
transmission. This mass was located on the frame at the center of 
gravity of the engine. The motor mounts were assumed to be made 
of rubber and the equivalent viscous damping for them was taken to 
be .05 times the critical damping for the engine mass-spring sys- 
tem, G, = 0.05 G, <rit = 0.1 K,M,. This value of damping corre- 
sponds to a decrease in amplitude of free vibrations of 30 per cent 
per cycle. (4) The system was assumed to be linear and the tire 
damping was neglected. 

Since, it was expected that any effect of the engine would be 
strongest near the natural frequency of the engine on its mounts, the 
response of the car to travel over a continuous wavy surface, such 
as a washboard road, was studied. Observations of the double am- 
plitude of the midpoint of the chassis were made for different engine 
mount stiffnesses and travel speeds. The midpoint of the chassis is 
the approximate location of the front seat of the automobile. 
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JOHN A, CARLSON 


An electrical analogy which was used to simulate the automobile 
is shown in Figure 1 also. (5) Other types of computers which could 
be used to solve the differential equations of motion for the me- 
chanical system are the differential analyzer (6,7) and digital com- 
puters. (8) Also, in this case, classical Laplace transform tech- 
niques could be used. 

The input functions for the electrical analogy were generated by 
an oscillator. These were sinusoidal voltages which were propor- 
tional to the vertical velocities of the tire-road contact points and 
were applied to the circuit at e, and e,. 


RESULTS OF THE ANALYTICAL STUDY 


The curves presented in Figure 2 show the variation of the frame 
midpoint double amplitude with engine mount compliance at three 
speeds. Additional measurements were made of the amplitudes of 
the front wheel motion and the relative motion between the engine 
and the frame. The wheel amplitudes remained practically constant 
at each speed, independent of the engine mount stiffness, and were 
as given in Table 1. Significant values of the double amplitude of the 
relative motion between the engine and the frame for a speed of 13.6 
miles per hour are .23 inch for an engine mount compliance, 1/K, 
.14 x 107* ft/lb. and 2.3 inches for 1/K, = .20 x 10°‘ ft/lb. 


TABLE 1 


FRONT WHEEL BOUNCE 


Speed Double Amplitude of Front Wheels 
(miles/hour) (inches) 
13.6 5 
20.5 0.46 
27.3 0.18 
CONCLUSIONS 


As would be expected, the action of the engine mass-spring 
system as a dynamic vibration absorber on the vehicle system was 
greatest at a particular excitation frequency, that arising from going 











{ LSON 


mobile 
1 Could 
e me- 
l com- 

tech- 


ted by 
ropor- 
ts and 


frame 
three 
ides of 
engine 
mstant 
| were 
of the 
f 13.6 
1/K, 


pring 
1 was 
going 











AUTOMOBILE ENGINE AS VIBRATION ABSORBER 55 


over the wavy surface at 13.6 miles per hour, as seen from Figure 
2. This frequency, 10 cycles per second, was approximately the 
natural frequency of wheel bounce. it is apparent from the Table 
that the front wheels were in resonance with the excitation received 
from the road. The conclusion reached here is that the engine 
mount may be chosen such that the shake of the frame caused by 
wheel bounce may be reduced or be absorbed by the engine. Static 
deflection of the engine on its mounts for minimum motion of the 
frame under this condition was 0.117 inch, 1/K, = .14 x 10°* ft/lb. 
Greater values than this lead to larger motions of the frame nad 
excessive vibration of the engine, because of a resonance condition. 

For higher speeds, higher frequencies of road excitation, the 
peaks caused by engine resonance are much lower and occur at 
greater engine mount stiffnesses. The large change in amplitude 
with speed indicates the sharpness of the wheel bounce resonance 
peak which makes the use of a dynamic absorber more feasible. At 
some speed between 13.6 and 20.5 miles per hour the peak of the 
corresponding curve will occur at 1/K, .14x 10 * ft/lb. If this 
peak is too large a slightly lower engine static deflection may be 
chosen as a compromise. Lower frequencies of excitation below the 
wheel resonant frequency would result in lower amplitudes until the 
natural frequency of the sprung mass on its suspension was ap- 
proached. At this frequency the engine is effectively rigidly at- 
tached to the frame if its static deflectionis as indicated above. The 
fact that this motion is highly damped makes the use of a dynamic 
vibration absorber for this low frequency impractical. Besides, 
the required static deflection of the engine would be excessive. 

This study sets an upper bound on the vertical static deflection 
of the engine to avoid excessive vibrations of the engine and of the 
frame for road excitation. This limit may have to be compromised 
by considering the effectiveness of the engine mount in isolating 
other engine vibrations from the frame. 

The technique of analysis shown by the above example yields a 
knowledge of the performance of the system for a range of designs 
so that the most satisfactory may be chosen. It should be possible 
to use the engine mass to absorb other vibrations of the vehicle as 
well as wheel bounce. Probably the most significant modes of 
engine motion which should also be considered are pitch and roll. 

Using this technique the possibility of employing the engine 
inertia and mountings to absorb a torsional mode of vibration of the 
frame and body could be explored. At the same time the isolation of 
torsional vibrations of the engine might be studied. Reference 3 
Shows the results of using the engine to absorb a bending mode of 
vibration of the frame and body. 


JOHN A. CARLSON 
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Engineering Simplifications Versus Mathematical Rigor 
for a Problem in Conduction Heat Transfer 


By Thomas D. Shea, Jr. 
Research Laboratories 
General Motors Corporation, Detroit, Michigan * 


ABSTRACT 


In many engineering problems where mathematical representa- 
tion can be derived in the form of equations and boundary conditions, 
developing a scheme for the solution of such a system is normally 
an arduous task even withavailable large computers. Such an analy- 
sis is often inconsistent with the time available to obtain a solution 
and the accuracy of knowledge concerning the problem parameters. 

For these reasons the engineer commences to introduce a series 
of simplifying assumptions to reduce the complexity and time re- 
quirements of the problem. The justification of these assumptions 
are founded principally on familiarity with the physical problem. 

In this paper the author deals with a problem in conduction heat 
transfer; namely, the determination of the temperature distribution 
ina radially symmetric flat disk in which the edge is held at con- 
stant temperature and cooling takes place on the sides by a medium 
at constant temperature. By means of a numerical example it is 
shown that the agreement is good between the simplified approach 
and the exact mathematical formulation. 


In many engineering problems where exact mathematical repre- 
sentation can be derived in the form of equations and boundary con- 
ditions, it is frequently an arduous or impossible task to solve for 
the variable of interest explicitly, thus to obtain a solution one must 
resort to numerical methods. Again, if one is familiar with the 
physical problem a numerical scheme can normally be devised that 
is more direct and simple which adequately describes the problem. 
To demonstrate this approach, a problem from conduction heat 
transfer will be investigated. 

Consider a radially symmetric disk, that is also symmetric in 
the plane of the disk, as shown in Figure 1. Because of symmetry, 
only one quadrant of the r-w plane is shown, The half-width, which 
is a function of radius alone, is given by w(r) as indicated, On the 
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rim r = b, let the temperature be constant and denoted by T), and let 
the disk be cooled on the sides by a medium at constant temperature 
Tc. Tofind the temperature distribution in steady state, the prob- 
lem will first be solved neglecting the temperature variation per- 
pendicular to the plane of the disk. 


THE ONE DIMENSIONAL SOLUTION 


To derive the governing differential equation consider a generic 
element of the disk in Figure 1. For steady state conditions the heat 


ae 


flowing into this element must equal the heat leaving.* If it is as- 


sumed that T, , gy > Ty > Tc, then the heat flows will be positively 
directed as shown in Figure 2, Because of radial symmetry the net 
heat flow in the 6 direction will be zero, The energy balance for 
this element is 


(1) Qr + dr = Qr + 2Q 


Using Fourier’s heat law which relates the heat flow tothe tempera- 
ture gradient equation (1) can be written 


(2) 2k-(redr).w (redr) (25) do - 2krw(r) (27) ao 


= 2h(T-T,) [ “°" rdsde 


, 


where k is the thermal conductivity of the element and h is the 
heat transfer coefficient. The line integral appearing above can be 
transformed to an integration with respect to r as follows 
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by their Taylor series expansion and the indicated products formed, 

| then in the limit as dr approaches zero, equation (2) reduces to 
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| where the second order and higher differentials in r has been neg- 
| lected, 

For a few particular functional relations between width and ra- 

dius the differential equation can be solved for the temperature 

| variation. One such relation occurs when the disk is flat. Equation 

| (3) can be reduced to the modified Bessel equation 
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where N = (h/kw]*/? and v is related to the temperature by the) yhere t 
transformation 





' 
ve=T-T. 
With the above relation the boundary conditions become =, 
, alt 
' 
(5) vy > Tp - Ig = Vb 
r=b (8 
and 
where 
dv _ aT} i the sul 
(6) drir = drirs0”° a sim: 
heat lc 
where the last boundary condition expresses the fact that for this 
radially symmetric problem, the heat flow at r = 0 is zero. ( 
The general solution of the differential equation (4) is | 
v = C,1,(Nr) + C,K, (Nr) where 
the b 
Evaluating the arbitrary constants by (5) and (6) yields heat | 
; anar 
(7) vey 2 odin 
” 1, (Nb) | taine 
rim t 


where the expression involves hyperbolic Bessel functions of zero be su 
order, Hence equation (7) represents the solution for'a flat disk T, is 
subject to a temperature T), on the rim and cooled on its sides bya the c 
medium at temperature T<. If the parameters of the problems are more 
specified, the temperature distribution can be obtained from equa- 
tion (7) by means of tables containing these computed functions,’ 
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THE ONE DIMENSIONAL FINITE DIFFERENCE FORMULATION 
] 
To develop a finite difference scheme that is not as severely re- prot 
stricted (one that can be used on a disk of arbitrary radial contour) } asa 
let equation (1) be written in the form accc 
Or’r+dr 2k(r + dr) w(r + dr) d@ : ey 
be ¢ 


For additional information see G. N. Watson “A Treatise on the 


Theory of Bessel Functions,” 2nd Ed., Macmillan Co., New York 
1944 
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where the derivative of temperature can be replaced by 


oT) = T, + dr af T, + € 
OP's s ay dr 


( 





in the limit, as dr approaches zero e¢ will approach zero. Introducing 
a finite number of stations, neglecting higher order effects, yields 


Qn + 2Qns 


4mr,+, Wn , 





(8) T,+, = Tn + (To +, 7° Fn) 


where the quantities Q,, and Q,.,. are defined as heat flows through 
the surface r =r, and the annulus surface between r,,, andr,. Ih 
a similar manner as employed in equation (2) for determining the 


heat loss to the cooling medium, Q.,, in equation (8) can be written 


5 + Fe 


rn + 
(9) Qsn = 2ah(T,, = To) ( 2 


) Cy 


where C,, is the arc length between station n+1 and n. In view of 
the boundary condition of equation (6), corresponding to n=0 the 
heat entering Q, at this point of radial symmetry is zero. Now if 
an arbitrary temperature is assigned to T,, by means of (8) and (9) 
the temperature and the heat entering the next station can be ob- 
tained. Repeated application of these equations ultimately led to a 
rim temperature Tr which, in general, is not equal to T),. As might 
be suspected from the exact solution for the flat disk in equation (7), 
T, is related in a linear fashion to Tr so that in order to arrive at 
the correct rim temperature the process will have to be repeated no 
more than twice. 


ENGINEERING IMPLICATIONS OF THE FINITE 
DIFFERENCE SOLUTION 


In many engineering applications the present formulation of this 
problem is unrealistic because the cooling medium cannot be treated 
as a reservoir at constant temperature, indeed, often the cooling is 
accomplished by a flow of gas. Hence, if it is assumed that a one 
dimensional temperature distribution adequately describes the disk 
and cooling gas two simultaneous differential equations will have to 
be satisfied. In the forementioned finite difference scheme this can 
be accomplished by the formula 


T = T, 
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TEMPERATURE DISTRIBUTION ONE DIMENSIONAL EXACT 
AND FINITE DIFFERENCE SOLUTION 
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where Wis the weight flow and C, is the specific heat at constant 
pressure of the cooling gas. . 

In order to obtain an idea of the correlation between the exact 
one dimensicnal temperature distribution given by equation (7) and 
the solution from the finite difference method, it is necessary to 
specify the parameters in the problem. In later numerical evalu- 
ation these parameters will also be used; for an assumed flat disk 
b = 4 inches, w = 1/2 inch, T,, = 1400°F, T, = 350° F, h = 40.5 
BTU/(FT* HR °F) and k = 12 BTU/(FT °F HR). 

The temperature as a function of radius for these two cases is 
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indicated in Figure 3. For the finite difference solution the stations 
were placed at quarter inch intervals; even with this coarse sub- 
division the agreement is quite good. 


THE EXACT TWO DIMENSIONAL SOLUTION 


For the special case of a flat disk the exact two dimensional so- 
lution that accounts for temperature variation perpendicular to the 
plane of the disk as well as the radial variation is known. The gov- 
erning partial differential equation and boundary conditions for this 
formulation will be stated and the method of solution will be indi- 
cated. For steady state conditions, the Laplacian in cylindrical co- 
ordinates is 


where again v is defined v = T-T,,. At the rim and origin the 
boundary conditions 


“lr=b- 


must be satisfied. Since the problem is symmetric with respect to 
the plane of the disk the following boundary condition can be em- 
ployed 


dvi 


Suin ee" 


On the cooling surface z = ¢ the variable v must satisfy the differ- 
ential equation 


where h is the heat transfer coefficient between the disk and the 
cooling medium and k is the thermal conductivity of the disk. A 
unique solution that satisfies the Laplacian in the disk and the 
boundary condition on the surface is 


Z 
Cos(y,, y? 
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where c = h/(kf/) and the values of y are the positive roots of ty 
equation® 


yn tanyn =c 


COMPARISON OF THE ONE AND TWO DIMENSIONAL 
SOLUTIONS 





' 
The temperature variation with radius on the plane of syuail 
z = 0 and at the face z = where heat is being extracted from the} 
disk are shown by the solid lines of Figure 4. Also shown in Figur} 
4 by the dotted line is the one dimensional, finite difference solution} 
It can be seen from this plot that neglecting the temperature vari-| 
ation in this problem, where the width to radius ratio equal o 
quarter, the deviation is small, 
The prominent reasons for the successful determination of a} 
exact two dimensional solution where the flat disk geometry, a heat} 
absorbing medium of constant temperature, and a constant hed! 
transfer coefficient. 
In engineering application some or all these factors migh} 
change significantly with radius. These, of course, can be readily 
accounted for in the one dimensional finite difference scheme. Sine 
only flat disks have previously been compared, the question arises 
is there a method by which the one dimension finite difference 
scheme can be checked when the disk is of variable thigkness ? 


NUMERICAL SOLUTION OF THE TWO DIMENSIONAL PROBLEM | 


For this case the one dimensional finite difference formulation 
can be compared to the two dimensional solution obtained by a re- 
laxation method. This scheme is simply a numerical method fo 
obtaining a two dimensional solution, To demonstrate the procedur 
and to indicate its sensitivity, the flat disk problem will be solvet 
by this method. As stated earlier the temperature distribution must 
satisfy Laplace’s equation which can be written 


oy og? at 2 op OT), at 
V r or dr 0z* 
2 : , 
The derivation of the two dimensional! solution is carried out! 
amanner similar to that used by Carslaw and Jaeger, Conduction 
Heat in Solids,” Oxford Univ. Pre , 1948, pp. 183-90. The solut 
iven here is somewhat simpler because there i ymmetry with re 


pect to the plane of the disk, 
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Now let (n,m) be an interior point of the disk as indicated in Fig- 
ure 5. Also iet the mesh size be selected such that Ar=Az. The 
relaxation procedure is initiated by assuming a temperature dis- 
tribution throughout the disk so that, in general, V°T # 0 at a ge- 
neric point. If the reSidue R, |, is defined as 


r 2 2 
Ram = (42)? 9° Tam 


n,n 





then in view of equation (10), the residue can be written as 
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(m) 


n-l,m 








y 


where rn+,/2 is the average of rn+,; and ryn,andr,-;/2 is defined in 
a similar manner. The boundary points on the z = 0 plane can be 
treated as an interior point, because of symmetry with respect to 
this plane, if the equation 


Pn+i/2 P Pn-1i/2 . — =: 
Ry om ” _ To m+ ad rn Tym-1 * - Th om - 4T,, 
is used instead of (11). On the axis of radial symmetry the boundary 
condition 


Ta,0 "7 Ty J 


must be satisfied. If the maximum value of n is N corresponding to 
the surface z = ¢, the equation 
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TEMPERATURE DISTRIBUTION ONE AND TWO 
DIMENSIONAL NUMERICAL SOLUTIONS 
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represents the boundary condition on the cooled surface. Again, if 
M is the maximum value of m, then 


on the rim of the disk.” For a fixed grid size the accuracy of the re- 
laxation procedure depends upon reducing the residue at each in- 
terior point to some epsilon sufficiently close to zero. 

The variation of temperature with radius for the one dimensional 
finite difference and the two dimensional relaxation method are 
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plotted in Figure 6. The temperature distribution of the plane ¢ 
symmetry z = 0 and plane of cooling z = ¢ are indicated by the solig 
lines. The mesh size, A z, used was one-quarter of an inch. Again, 
it is seen that the one dimensional solution is sufficiently close ty 
the numerical two dimensional solution for many engineering ap. 
plications. The discrepancy between the two dimensional exact so. 
lution of Figure 4 and the solution shown in Figure 6 obtained by the 
relaxation is due to the difficulty of satisfying all of the boundary 
conditions simultaneously on a grid this coarse. 


CONCLUSION 


In the above discussion, the heat transfer problem of determin- 
ing the temperature distribution for the disk was approached from 
four points of view. In addition, a comparison was made of thes: 
approaches by means of a numerical example, Although there is 
some latitude in the selection of the method, the approach chosen 
must be consistent with the application and the intended use of the 
resulting information, In many engineering applications the simpli- 
fied one dimensional! finite difference method gives adequate repre- 
sentation of the temperature distribution. Additional accuracy is 
usually unwarranted because of greater discrepancies between the 
physical problem and the mathematical model, or because of the un- 
certainty in determining accurately the heat transfer coefficient 
when gas at high Reynolds numbers is the cooling agent. 
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Annuity Interest 


By Wright Collier 
Consulting Engineer 


Introduction. 


The engineer’s concern with the problem of calculating annuity 
rates is twofold. In his professional capacity he encounters it in 
such applications of engineering economics as setting up amortiza- 
tion, depreciation, and equipment replacement schedules and in the 
more or less continuous routine of comparing proposals. As a pri- 
vate citizen, his interest in the problem has been heightened by the 
great expansion, in recent years, of the deferred or distributed pay- 
ment idea, 

These notes present an iteration method forthe rapid calculation 
of annuity interest rates with a relative error 0.01. The method 
combines two familiar techniques: guess-work and log-log slide- 
rule calculations. 


Notation and Equations. 








For an ordinary annuity: 


n = number of payments outstanding. 
m = number of payments per year. 

r = interest rate per payment period. 
x = atrial value of r. 


i = effective annual interest rate. 

R = annuity rental = payment to be made at the end of each of the 
n payment periods. 

A = present value at interest r, rental R and with n payments out- 
standing. 

a% = present value at interest r, with n payments outstanding and 
rental = $1. 

a, = value of a, for r = x. 


relative error in calculated effective rate, i. 
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isis (4) 


Application. 


The method consists of solving equation (3) for each of a series 
of values of x. The trials are continued until a value of x is found 
for which ax=a, This value of x=r. Equation (4) is then solved 
to obtain the effective annual rate, i. 

There are some simple aids to the formulation of a rapidly con- 
verging sequence of guesses at the value of r: 


1. If a current local rate, or range of rates, for the type of 
transaction under consideration is known beforehand the work of 
computation may be considerably abridged. Lacking this informa- 
tion, the number 0,01 has several features that recommend it as a 
first guess at the value of r. It is reasonably central with respect 
to the rates encountered in practice, and a close approximation te 
the expression 1.01" can be evaluated directly on the rule, fora 
wide range of exponents, without setting the slide. Also, 9.01 makes 
a convenient denominator for equation (3). This value will be used 
first in the examples that follow. 


2. In any cycle of the solution, if a > a,, the value of x is too 
small. Conversely, if a, is too small, x is too large. 


Example 1. 


The price of a record player is $300. The terms of Sale are $50 
down and payments of $15 at the end of each month thereafter, for 
eighteen months. What rate of interest is being charged? 


A = $250 ao 250/15 
R $15 $16.67 
n = 18 
m 12 
x ies 1 - (1 + x)" ¥ 
0.01 0.836 0.164 16.40 


This is fairly close for a first approximation. Try 


0.008 0.866 0.134 16,7+ 
Try 
0.009 0.851 0.149 16.6- 


The value of r lies between the last two values of x and is somewhal 
closer to 0.008. Take r = 0.0083. 
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0.0083‘ - 1 = 10.4 per cent 


A parallel calculation of the entire problem, using 7-place loga- 
rithms, gives i 10.322 per cent. Hence, 


400 - 322 


10322 0.0075 


Example 2. 


An automobile dealer offers a late model of a certain make car 
at $1595. Terms of sale are $150 down and $9.97 at the end of each 
week after date of sale for a period of three years. What is the an- 
nual rate of interest on this contract? 


A $1595 - $150 $1445 R = $9.97 
n = 156 a = 1445/8.87 
m 52 = $144.94 


The solution of the preceding example was given in detailed tab- 
war form. Actually, the solution of these problems requires no 
writing — evaluation of the expression 1 - (1 + x)" is an easy 
mental operation; the rest of the work is mechanical computation. 
In this, and in the final example, only values of x and the corre- 
sponding values cf a, will be recorded. 

Xx a. 
0.01 78.8 
This result is so far off that a drastic reduction of x is in order. 
Try 


0.001 144 
Try 

0.0008 146- 
Try 

0.0009 145+ 


r lies between 0.0009 and 0.001. Take r = 0.0094. Then 
i 0.0094 - 1 4.97 per cent 


Solution with logarithms gives i 





5.005 per cent. Therefore 














Example 3. 


A debt of $250,000 is to be retired by payments of $1,000 at the 
end of each month after its inception. Payments are to continue for 
99 years. What is the interest rate? 


A = $250,000 R = $1,000 
n = 1188 a, = $250 
m 12 
x ax 
For x = 0.01,(1 + x)" 0.0000074 and 1 - (1 + x)" = 0.9999926, 
If it is assumed that this numerator will remain fairly close to 1,0 
for the remaining cycles, then, for x 0.004, a, should be approx- 


imately equal to 250, 


x _ay 
0.004 248- 
Try 
0.0039 254- 


These are bracketing values and r is closer to 0,004, 
Use r = 0.000397. Then 


i = 1.00397 - 1 4.87 per cent 


Recalculation with 7-place logs gives i 4.886 per cent. So, 
16 
€ 4886 0.003+ . 


It is worthy of note that, as n becomes large, the interest rate ofan 
ordinary annuity converges to that of a simple perpetuity; that is, 


LirT 


r R/A. For this example, the limit has the value 1,000 


— @ 


250,000 = 0.004. This would have been the logical first choice for 
the value of x, 


Conclusion, 


The discussion here has been restricted to the basic (ordinary) 
annuity. More complicated annuities still have more complicated 
solutions. It is not difficult to extend the principles of this method 
to calculation of the interest rate of other types of annuities. 


Reference. 





Harper, Floyd S. “Mathematics of Finance” — first ed. 
International Text Book Co., 1946. 
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The Problem of the Hampered Input 


By Federico Gabriel 

Systems Analysis Laboratories 
Hughes Aircraft Company 
Culver City, California 


ABSTRACT: 


Given a set of interrelated inputs (with storage deposits) and 
outputs, and a situation in which an input is forced to supply below 
its maximum for a certain interval of time, the following problem is 
studied: What is the optimal rate of supply for every other input ? 
The general solution ot any problem of this type is presented in 
terms of a special matrix representation. 


The problem to be studied originated from the following situ- 
ation. A factory in Argentina had four lines, each with the corre- 
sponding storage deposit at the end, producing parts that entered in 
the production of two subsystems; the subsystems having line No. 3 
incommon. Line No. 2 suffered damages that, according to the 
maintenance and repair department, would bring down the production 
rate of that line to about a third of maximal output for a period of 
twelve to thirteen weeks, until spare equipment arrived from Europe 
and was installed. The storage deposit for that line, together with 
the concurrent output of the damaged line, could not maintain the 
feed required for maximal output of the subsystem. The manage- 
ment was faced with the question: What is the optimal rate of pro- 
duction for each line ? 

The problem is one of supply and distribution, and the lines and 
Storage deposits may be considered as inputs to the subsystems 
production. Therefore, the problem is of inputs and outputs, when 
one of the inputs is supplying below its maximum. 

The accompanying diagram illustrates the geometry of the 
Situation encountered: 


Notation: 
I; = Input No. j (j = 1, 2, 3, 4) 


0; = Output No. j (j = 1, 2, 3, 4) 


I; = Rate of production of I; during maximal operation. 











O Rate of production of O during maximal operation. 
D Storage deposit No. j (j 1, 2, 3, 4) 

F Full storage deposit 

O Empty storage deposit 

T Duration of the period of reduced production 

D 


Initial boundary condition = quantity of output of I; in storage at 

the beginning of T. ® 

D Final boundary condition = quantity of output of I; in storage at 
the end of T. 

[' Rate of production of I; while damaged (during T). I' < I 

O' - Rate of production of O during T. 


Intersection of 1 if they intersect (i.e., if I; feeds Ox) 
I; with O O if they do not intersect 


The weights of each input towards each output are as shown in 
the following table (these weights measure the number of unit inputs 
per unit output): 

At instant t, input I, suffers a disruption in its mechanism that, 
for an estimated period of T units of time, will decrease its oper- 


y 
ation to a fraction - of its maximal performance. What adjustments 
Ss 


should be made in the other three inputs inorder to maintain optimal 
operation during the period of decreased production? 
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It may continue 
working at maximal production. In fact, O, output will remain maxi- 


The solution of the problem is immediate for |: 
mal, since the consequences of the emergency on I, , if any, will be 
to decrease the fraction of its production absorbed by O,. Hence O, 
will keep receiving its full share of I, towards maximal production. 

Whether or not the emergency will affect I, at all depends on 
whether D, - D, is or is not smaller than the difference between the 
maximal output of I, during a normal period of duration T and that 
part of it required by O, and O, combined during the emergency 
period T. 

In order to use the cushioning properties of the deposits to the 
utmost we must have the following boundary conditions: 





| f 
Pe 2 D 
I, | ad. ' Table of 
I. | D, O Boundary 
| . 
[, | pi k Conditions 
[, | D. | D, = Dy 


Since I, and I, may have to be regulated, we need to know the 
effects of such regulation ppon the efficiency of the system, in order 
to regulate in optimal manner. For that we need the functional 
relation between production efficiency and input quantity. Assuming 
anon-linear and monotone relation, it can be demonstrated that the 
Optimal regulation of input consists in a constant rate during T, 
rather than in any monotone or oscillatory function of time.’ 
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Our next step is to set up a material balance for the period T, 

From the table of weights we see that one unit of O, is in corre- 
spondence with P,, units of I,, with P,, units of I, and with P,, units 
of I,. An analogous set of relations holds with respect toO,. We 
may represent these correspondences as follows: 


O, < > Pi, I. P., I.. Ps, I, 


O, * > Po ly Po ly 


Considering the boundary conditions follows that the contribution 
—_ Vv 
of L, towards O, during T will be D, + : I, T. Applying (1), the 
_ . V ‘ 
production of O, during T will amount to P,, (D,+- I, T). 
1 2 ha Sex 


On the basis of (1) and of the principle of constant rate of input 
mentioned in a previous paragraph we are now able to obtain the 


0 
optimal production rates to secure of I,andI,. In effect, I,< 2 


31 
and using the indicated boundary conditions we see that, during T, 
we have 
i p 
(F - D,) +5” (D, + 


“ . I, T) units as the total for I, . 
31 


On the other hand, optimal production of 0, during T requires 


1 units of I,. Therefore, since I, must supply both O, and . 
O.T I Pp i : 


- +(F -D)+ p- (D+ : I, T) units as the total for I, during T. 
" Dividing the last two expressions by T we obtain the production 
rates to establish for I, and F 

With the use of matrices these analytical results may be put in 
a form suitable for immediate extension to any number of lines and 
different interconnections between inputs and outputs. 

As a first step the above expressions are rewritten in tabular 
form, with the understanding that the sum of all the elements of row 
no. j gives the rate of production contributed by I; during T: 

Next the notation for the table of weights is generalized: Pix 
shall relate the number of units of I; that correspond to one unit of 
Ox . 


If the damaged input is I, and it feeds O, , then I! I, and O} = 


V 
Ss 


a tr + 2. 
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PROBLEM OI THE HAMPERED INPUT 


THE 





Contribution Cont ribution Contribution 
to Storage to O, to O, 
y 
- F - D, P,, (D2. +- I, T) 0, 
ee s : 
T Pia © P,. 
—+— + - 
, D, ! ; : 
I, - + ;- O Per Unit Time 
: T ao 
——_+— + + 
F -D, P, (Di +* I, T) 
~ Ss 
I 
; r Pon 2 o 
I O O | ©, 
4 Pp , 
Pa 42 


Introduce the more general notation in the above table and de- 
compose into matrices as follows: 


Cy C., C;, C,, 
/ 
| D)-D A a ee 
! ! l 2 2 2 
ee Oe Me 
C., C., C;, C4, 
Pi, P.. Ps, Pa. 


where the index j in the first element of the first matrix takes the 
value of the number of the column that it multiplies on the second 
matrix (the purpose of this abbreviation is to avoid writing a sum of 
four products of matrices 4 x 3 by 3 x 4). 

The above product yields four sums, one per each input. The 
second matrix takes care of the interconnection geometry and of 
localizing the damaged input; the first one accounts for the “sinks” 
of the inputs (borrowing an expression from the theory of vector 
fields). Notice that D, is a “negative sink” or “source,” since D, 

0 and D, 4 0. 

The mode of representation (2) is applicable to any number of 
inputs and outputs, and to any localization of damaged inputs. Any 
One specific problem is individualized by four items; The inter- 
connection numbers, the table of weights, the table of boundary 
conditions, the localization of the damaged inputs in the matrix on 
the right in (2). Matrix representation (2) illustrates how these four 








items are combined to obtain the solution to any problem of this 
type, in an almost automatic way. 

In case that, before the end of the period T, conditions change 
in such a way as to indicate that the emergency situation is to be 
prolonged beyond t, + T, the method presented in the paper referre¢ 
to in the second footnote can be effectively employed to handle the 
new situation. 

Naturally, the determination of the storages Dj; that are optimal 
may be considered, in itself, an operations research problem during 
the design stage; as also will be the optimal method for returning 
the reservoirs from D;, to D; after the emergency period (this 
optimal method is a function of the probability of damage at I; and of 
the relation between total production efficiency and I;). 

Should one desire to extend the above results to more complex 
types of interconnections (where outputs become inputs of further 
outputs, etc.) the Kron approach to complex systems may be enm- 
ployed. ° 


Federico Gabriel 


APPENDIX 


Denote by C(R) the unit cost corresponding to a rate of produc- 
tion R (for example, to a daily rate of production). Assume that the 
functional relation between those two variables is of atype repre- 
sentable by the curve in Figure 1. 

Determining the optimal regulation of input consists in estab- 
lishing R as a function of time t such that 


>, -T 
/ C(R) R(t) dt is a minimum and V J R(t) dt, where V 
O O 
is the total volume of production during T. 
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Figure | 


This is an isoperimetric problem in the calculus of variations, 
with variable end points. To solve it with the techniques of the 
variational calculus, continuity of R as a function of t must be as- 
sumed. The desired functional relation, however, is accessible to 
more modest mathematical tools; and the proper mathematical 
formulation will, collaterally, free the study from the idealization of 
continuity. 
dR 
(the change in monotonicity for values of R beyond the optimal may 
rev be neglected, since they have no bearing on the problem at hand), 
and because any rate of production that might be used would be 
maintained for a finite interval of time (i.e., the rate of production 
is not changed continuously, in the mathematical sense of the word 
“continuous”), the problem may be expressed in terms of finite 

John increments as follows: 
tatir Comparison is to be made, from the standpoint of total cost 
pan during time At, of the two following courses of action or production 
Schedules: see Figure 2. 


In fact, because | is a monotone decreasing function of R 


A - Constant rate of production R, during At. 


MET B - Rate of production R. +A R during “ , 


A 
and rate of production R,. - AR during , 
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A 


At ‘ , 
Here — is perfectly arbitrary,thus making the solution adapta- 


ble to the different possible circumstances in practice. On the other 
hand, the total volume of production during At must evidently be 
equal for both courses of action. 
ao dC | 
Applying the decreasing monotonicity of | to the given 
formulation of the problem the following relation holds: 


R . + A 2 cA R « (Al) 
c c 


To simplify the algebraic manipulations, the following abbrevi- 
ations will be introduced: 


C; total cost, during At, using production schedule A. 


i 


S. = total cost, during At, using production schedule B. 


Cr will be denoted by C, . 


c 


Cr - AR Will be denoted by C,. 


Cr 4 AR Will be denoted by C, . 
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Performing the necessary operations: 


—_— At ra 
€,-E.=[C,(R. - AR)- CR] F -[ CR - CR. + AR)] F = 
: At 
=[ C,(R. - OR) - C,(Re - AR) +C,(R. - OR) - C\Re | + 
At 
-[ C,R. - C,R. +C,R. - C,(R. + AR) | 5 te 
= [ (C,-C\)R.-(C,-C,) AR-C, AR-(C,-C,)R.+C, AR] & = 
, At 
“ [ {(C,-C,)-(C,-C,) } R, - {(C,-C,)-(C,-C,)} AR | — 
At 
=[ (C,- C,)- (C,- C,)] [Re - AR]. (A-2) 


Because R. > AR, inequality (A-1), implies that expression 
(A-2) is larger than zero. In consequence, production schedule A is 
optimal. 
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